Stability of spinning ring solitons of the 
cubic-quintic nonlinear Schrodinger equation 



Isaac Towers 1 Alexander V. Buryak x , Rowland A. Sammut 1 , 
Boris A. Malomed 2 , Lucian-Cornel Crasovan 3 and 
Dumitru Mihalache 3 

1 School of Mathematics and Statistics, 
University of New South Wales at ADFA, Canberra, ACT 2600, Australia 

3 Department of Interdisciplinary Sciences, Faculty of 
Engineering, Tel Aviv University, Tel Aviv 69978, Israel 

3 Department of Theoretical Physics, Institute of Atomic 
Physics, P.O. Box MG-6, Bucharest, Romania 



Abstract 

We investigate stability of (2+l)-dimensional ring solitons of the nonlinear Schrodinger 
equation with focusing cubic and defocusing quintic nonlinearities. Computing eigen- 
values of the linearised equation, we show that rings with spin (topological charge) 
s = 1 and s = 2 are linearly stable, provided that they are very broad. The stabil- 
ity regions occupy, respectively, 9% and 8% of the corresponding existence regions. 
These results finally resolve a controversial stability issue for this class of models. 
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1 Introduction 



Recently, much interest has been focused on spinning optical solitons, i.e., 
those carrying topological charge, in both (2+l)D [1-5] and (3+l)D [7-10] 
geometries. A spinning soliton has an embedded phase dislocation and carries 
intrinsic angular momentum. The integer number of phase rotations around 
the dislocation is the soliton's topological charge or "spin" . 

Broadly speaking, spinning solitons can be divided into two classes: (i) dark, 
i.e., vortices produced by a phase dislocation which is embedded in an infi- 
nite background; and (ii) bright, with the vortex core embedded in a bright 
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(localised) multidimensional soliton proper; the amplitude of which vanishes 
at infinity. In this Letter we consider bright spinning solitons in the (2+l)D 
geometry. In terms of nonlinear optics, (2+l)D solitons may be naturally re- 
alised as spatial solitons in the form of cylindrical beams in a bulk medium, or, 
alternatively, as spatiotemporal solitons in the form of fully localised "light 
bullets" propagating in a planar waveguide (film). Due to the presence of 
the vorticity, the soliton's cross section has an annular shape. We shall refer 
to them in what follows below as localised vortex solitons (LVSs). As for the 
(3+l)D solitons, they are spatiotemporal "bullets" propagating in bulk media. 

Models that may give rise to stable multidimensional, (2+l)D or (3+l)D, 
solitons must necessarily have nonlinearity which prevents dynamical collapse. 
Well-known examples of collapse-free nonlinearities are x^ (second-harmonic- 
generating), saturable, and cubic-quintic (CQ); the cubic and quintic compo- 
nents being, respectively, self-focusing and self-defocusing. All these nonlin- 
earities occur in various optical media. 

However, unlike ground-state (zero-spin) (2+l)D bright solitons, for the spin- 
ning ones the absence of collapse does not guarantee dynamical stability. In 
fact, a LVS tends to be strongly destabilised by azimuthal perturbations which 
break it up into several separating zero-spin bright solitons (the latter are 
stable). In (2+l)D models with x^ an d saturable nonlinearities, numerical 
simulations had revealed strong azimuthal instability [1,2], which was later 
observed experimentally in a x^ medium [3]. A similar instability of (3+l)D 
rings resulting from the x^ nonlinearity has been found in Ref. [10]: a 3D 
LVS (in fact, it is a torus) breaks up due to an azimuthal perturbation, and 
the resulting zero-spin solitons fly off in directions tangential to the ring. 

In terms of the LVS stability, more promising are models with competing non- 
linearities. The first step in this direction was the study of (2+l)D rings in the 
CQ nonlinear Schrodinger equation. Simulations reported by Quiroga-Teixeiro 
and Michinel [4] showed that broad rings found (by means of the variational 
approximation and direct numerical methods) in that model not only did 
not demonstrate growth of small perturbations, but also survived collisions, 
thus appearing fairly stable. However, a similar analysis for the (3+l)D CQ 
model [9] has demonstrated that narrow LVSs were broken up quickly, while 
broad rings were destabilised much slower by azimuthal perturbations. How- 
ever, eventually all the LVSs for which a definite numerical result could be 
obtained were unstable. This suggests re-checking the above-mentioned sta- 
bility of LVSs in the 2D model reported in [4] . Rerunning simulations for the 
same cases which were considered in [4], it was established [6] that, in fact, 
they are also subject to the weak instability against azimuthal perturbations, 
provided that the simulations are long enough. 

As 2D and 3D LVSs in the CQ model become very broad, it still remains 
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unknown whether the growth rate of their instability against azimuthal per- 
turbations gradually vanishes in the limit of infinitely broad rings, or if there is 
a clearly defined transition to truly stable solitons. The problem of discerning 
between very weak instability and true stability may be of little importance for 
applications, as experiments are always carried out in finite samples, which do 
not have enough room for the development of an instability if it is extremely 
weak. Nevertheless, the issue is of principal interest and is therefore worthy of 
consideration. 

Another model with competing nonlinearities which may be promising for the 
generation of stable rings combines the x^ an d self-defocusing cubic [x^] 
nonlinearities. It is necessary to say that no conventional nonlinear material 
with strong x^ nonlinearity directly satisfies the requirement of this model to 
have a negative x^ coefficient at both the fundamental- and second-harmonic 
frequencies (see below). However, two possibilities to create the necessary ef- 
fective x^ nonlinearity have been proposed: (i) by creating a layered medium 
in which layers providing for the x^ nonlinearity periodically alternate with 
others that account for the self-defocusing Kerr nonlinearity, and (ii) by en- 
gineering special x^ quasi-phase-matched gratings [12]. In the latter case, 
induced an d intrinsic x^ nonlinearities may be equal in strength, and 
the former one may be given either sign. 

Recently, LVSs were considered in this x^ ~ model [11]. Using direct 
simulations and linear stability analysis, it has been shown that narrow rings 
demonstrate typical breakup into zero-spin bright solitons initiated by the 
azimuthal instability, but very broad (flat-top) LVSs with the spin (topological 
charge) s — 1 and s = 2 were indeed found to be completely stable. In these two 
cases, the stability region is, respectively, ~ 8% and ~ 5% of the corresponding 
existence domain. 

In the cascading limit, corresponding to large wave- vector mismatch, the x^ — 
■^(3) model reduces to the CQ model, which suggests that there may be a chance 
to find completely stable LVSs in the CQ model too, which is a subject of the 
present work. We will employ the same techniques that were applied in Ref. 
[11] to the x^ model, i.e., rigourous computation of the stability eigen- 

values. This will overcome limitations of "phenomenological" analyses of the 
CQ model carried out in previous works, which relied on simulations of per- 
turbed solitons to directly test their stability, or simulations of the linearised 
equation to estimate the largest growth rate of instability. A shortcoming of 
both methods is that they can miss very weak instability. We will demonstrate 
that LVSs in the CQ model are truly stable if their width (or energy) exceeds a 
certain critical value. This result completes the investigation of the CQ model 
in the (2+l)D case, which was the subject of many above-mentioned works, 
and suggests a challenging question: can spinning "bullets" with sufficiently 
large energy be stable in the same model in the (3+l)D case? 
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Fig. 1. The stationary (2+l)D spinning-soliton solutions: (a) s = 0, (b) s = 1, and 
(c) s = 2. The values of k are indicated near the curves. 

2 The ring solitons 



Following the derivation of Ref. [8], but assuming two transverse spatial di- 
mensions, we arrive at the dimensionless CQ equation 

i— + V 2 { u + \u\ 2 u - \u\ A u = 0. (1) 
oz 



where u is the envelope of the electromagnetic wave propagating along the 
^-direction in the optical medium. In the case of the cylindrical beams in the 
bulk medium (i.e., spatial solitons, see above), the Laplacian in Eq. (1) is the 
diffraction operator acting on the transverse spatial coordinates x and y. 

In the alternative case of spatiotemporal solitons in the planar waveguide, y is 
replaced by a properly scaled temporal variable, r = t — z/Vq, where t is time, 
and Vq is the group velocity of the carrier wave. In the latter case, the part 
d 2 /dr 2 in the Laplacian accounts for the temporal dispersion (which must be 
anomalous, in order to have to right sign), rather than diffraction, and LVS 
will be, in unrescaled coordinates, a compressed (elliptic) ring moving in its 
own plane. 

Ring solitons are localised stationary solutions to Eq. (1) of the form 

u = U{r) exp(is9 + inz), (2) 



where r and 9 are polar coordinates in the (x, y) plane, k is a wave number 
shift (relative to the carrier wave), and the integer s is the spin. The amplitude 
U can be taken to be real, obeying an equation 

^L + ^JL^u. KU + u ,. v , =0 . (3) 



We solved Eqs. (3) by means of the relaxation technique. When k is small (a 
low-power regime), LVSs are narrow. The beam's amplitude at first increases 
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with k and then saturates, while the ring's width keeps increasing because of 
the self-defocusing effect of the cubic term. 

The wave number k parameterises a family of stationary solutions, examples 
of which are displayed in Fig. 1. The existence regions for solitons in the 2D 
and 3D versions of the CQ model are known [4,8]: 

< k < «Si « 0.18; < k < «22t « °- 15 > ( 4 ) 

and in both cases they practically do not depend on the soliton's spin [8]. 
The width and energy of LVS diverge as k — > K oSset . Note that a soliton 
solution to the ID version of Eq. (1) is known in an exact elementary form, 
the corresponding exact offset wave number being /^ffeit — 3/16 = 0.1875, so 
that the values (4) are close to it. 



3 Stability 

As a precursor to the full linear stability analysis of the ring solitons, we first 
consider the modulational stability of plane continuous- wave (cw) solutions to 
Eq. (1). This may give some insight into the stability of broad rings as they 
tend to these cw solutions in the limit k — > ^offset- The cw solutions with the 
propagation constant k are 

uq = a exp(inz), Oq = - (l ± y/\ — 4k) . (5) 

We take small perturbations to the plane wave of the form 

u\ = [a exp (ikx + iflz) + /Sexp (— ikx — iVt* z)} exp(inz), (6) 

where k is an arbitrary perturbation wave number, and Q is the stability eigen- 
value (the asterisk stands for the complex conjugation); the cw solutions being 
stable if Q is real for all real k. Substituting into Eq. (1) the perturbed solu- 
tion u = uq + ui and linearising around u , we find after some straightforward 
algebra that: 

n 2 = k 2 + k 4 ±k 2 VT^4^-4k 2 K. (7) 

It is easy to demonstrate that for the cw branch with the upper sign in Eq. (5) 
(i.e., with the larger amplitude), Q 2 > for all k , hence, this branch is always 
stable, while the other one is not. As it was said above, the LVS solutions 
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Fig. 2. Gray-scale contour plots illustrating the evolution of an unstable ring soliton 
with s = 1 and k = 0.05: (a) z = 0; (b) z = 310; (c) z = 330. 

of the CQ model that we are dealing with tend to the stable cw solution as 
ft —* ^offset, giving an initial indication that, as the rings broaden, they may 
become stable. 

This conclusion is suggested by simulations reported in Ref. [4], where suf- 
ficiently broad rings appeared to be stable, whereas narrower ones were def- 
initely unstable against azimuthal perturbations (see Fig. 2 ). On the other 
hand, as was mentioned above, more accurate (longer) simulations of the cases 
that were considered in that work show that weak instability still occurs. 
To perform the direct simulations we used the Crank-Nicholson scheme as a 
finite-difference approximation to propagation equations (1). The correspond- 
ing system of nonlinear equations was solved by means of the Picard iteration 
method ( see details in Ref. [13]). Typically, we chose equal transverse grid 
sizes Ax = Ay = 0.4 and the longitudinal grid size was Az = 0.02. Thus, 
direct simulations show a general trend of suppression of the azimuthal insta- 
bility with the increase of the size of LVSs. This may be sufficient to expect 
experimental observability of LVSs, but, following this analysis, one cannot 
predict if the instability may be completely eliminated, provided that the size 
of LVS exceeds some critical value. 

The stability issue can only be resolved in the precise sense by comprehensive 
analysis of eigenmodes of small perturbations around LVS. To this end, we 
add infinitesimal complex perturbations e(z, r, 6) to the stationary solutions 
of Eqs. (1) and (3) with the vorticity s, cf. Eq. (6), 

u — [U (r) + e(z, r, 9)} exp(is9 + inz). (8) 



A general perturbation e(z,r,8) may always be expanded into a series, with 
each term having its own vorticity J, so that a generic independent perturba- 
tion term is 

e = £+( r ) exp (i\ z + UQ) + £y(r) exp {-i\*z - J6) , (9) 
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where A is the (generally complex) LVS's instability eigenvalue. Substituting 
this into Eqs. (1) and linearising, we arrive at a non-self-adjoint eigenvalue 
problem: 



a 



+ 



D 



(10) 



—D -C. 



where £,= (#,£j),C± 



Lj + 2U 2 - 3U\ D = U 2 - 2U 4 , and 




— K. 



Instability is accounted for by eigenvalues with ImA ^ (the present system 
being Hamiltonian, eigenvalues appear in complex conjugate pairs or quadru- 
plets). The continuous spectrum of the eigenvalues consists of real intervals 
k < A < oo and — oo < A < —k. 

To analyse the eigenvalue problem (10), we replaced the differential operators 
by their fifth-order finite-difference approximations and solved the resulting 
algebraic eigenvalue problem numerically. We mostly used grids with 400 to 
800 points, but up to 1200 points were used in regions where a change of 
the stability occurs. To verify the precision of the numerical code, we also 
used another technique based on the relaxation method for solving two-point 
boundary- value problems. Although limited to finding real eigenvalues, the 
latter method admits a high degree of precision control without much of the 
computational overhead of other methods. For instance, it has been recently 
used to a great effect in finding a small stability window for higher-order soli- 
tons in a third-harmonic-generation model, which would have otherwise been 
overlooked [14] . Comparison between the spectral and relaxation methods has 
shown that the former one has good precision for the number of the grid points 
used: the numerical error in calculating the stability-boundary values of k is 
estimated to be 5n st ~ 10~ 5 for 1200 grid points. 

Results of the linear stability analysis for the fundamental (s = 1) and next- 
order (s = 2) LVSs are summarised in Figs. 3 and 4. We considered perturba- 
tions with J = 0, ±5, and have found that instability of the fundamental 
LVSs is not generated by perturbations with J > 3. The most persistent in- 
stability mode in the case s = 1 corresponds to J = ±2. Subsequent direct 
simulations demonstrate that this instability mode, if it takes place, initiates 
eventual breakup of the ring into several zero-spin solitons. In Ref. [2] an esti- 
mate for the number of filaments N resulting from the destruction of the LVS 
is given as N ~ 2s. Our results agree with this estimate: for LVSs with s — 1, 
dominant instability corresponds to J = ±2 and a singly charged LVS breaks 
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Fig. 3. Unstable eigenvalues for the ring solitons with s = 1. The spin J of the 
azimuthal perturbation is indicated next to each curve (the perturbations with 
J > 3 caused no instability, therefore they are not displayed). Only Im(A) is shown. 
Note that the dominant instability has J = 2, and it vanishes at k ~ 0.16, while 
the ring solitons exist up to k = K fj se t ~ 0.18. 
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Fig. 4. The same as in the previous figure but for s = 2. 

into two filaments, as can be seen in Fig. 2. (This agreement largely holds for 
the s = 2 case as well where for most of the unstable domain the dominant 
instability corresponds to J = ±4 and typically doubly charged LVSs break 
into four filaments.) 

In all the cases considered, we have found that there is a stability- change value 
K st , at which the largest instability eigenvalue ImA vanishes, and remains, 
along with all the other ones, exactly (up to the numerical accuracy) equal to 
zero in the stability window, tz st < k < K ofiset (recall K ff se t is the upper existence 
boundary of the LVS family); in other words, thin LVSs are unstable and broad 
ones are stable. The existence of the window is clearly illustrated by Fig. 3. 
For the fundamental LVSs (s = 1), the stability window occupies ~ 9% of the 
existence domain [0, ^offset] • 
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For the LVSs with s = 2, a similar situation takes place. The dominant in- 
stability in a larger part of their existence domain is generated by the pertur- 
bation with J = 4, but for broad LVSs it is overtaken by the J = 2 mode. 
The linear spectrum of the LVSs with s = 2 contains a stability window which 
occupies 8% of the existence region. 

It appears that stable LVSs cannot have the value of the spin larger than 2. In 
particular, the LVSs with s = 3 were found to demonstrate a persistent weak 
instability associated with the J = ±1 perturbation modes at all the values of 
k, and it seems very plausible that higher-order LVSs will continue to do so. 

In the work of Quiroga-Teixeiro and Michinel [4], the authors use variational 
techniques to make an estimate for the width of the stability window. Using 
the formulas derived in Ref. [4] the ratio of n st / K, oSsct for s = 1 and s = 2 LVSs 
is 0.803 and 0.838 respectively. In other words according to the variational 
analysis of [4] the stability window occupies approximately 20% and 16% of the 
existence domain [0, ^offset] for singly and doubly charged LVS. These estimates 
correctly predict the existence of stability domains for both s = 1 and s = 2 
LVSs and also the decrease in stability window size for the doubly charged 
LVSs relative to the singly charged. The factor of 2 difference between the 
estimates and the numerical results obtained in this work is quite reasonable 
considering the ansatz used in Ref. [4] becomes increasingly less accurate for 
large k i.e. where the stable LVSs exist. 

Lastly, it is relevant to mention that, in the work [6], a variational approach 
was used to study, in an analytical form, instability of LVSs in the (2+l)D 
and (3+l)D versions of the CQ model against a special perturbation mode 
in the form of an infinitesimal shift in the position of the vortex core relative 
to the broad soliton as a whole. In terms of the expansion (9), this mode 
corresponds to J = ±1. If 5H is a small variation of the value of the soliton's 
Hamiltonian generated by the infinitesimal off-centre shift of the vortex core, 
an instability condition which is well-known from the general soliton stability 
theory is 5 H < [15]. Assuming k close enough to ^offset, which is necessary to 
have broad solitons whose outer radius is much larger than the radius of the 
vortex core, it has been shown that a shift of the core leads to a decrease of an 
effective Hamiltonian of the interaction between the core and the outer rim of 
the soliton, which implies an instability. Because the interaction Hamiltonian 
is exponentially small in the case of the broad LVS, the instability is also 
expected to be exponentially weak. This instability mode was predicted for 
s = 1 and 2, but not for s > 3. It is noteworthy too that only for s = 1 the 
thus predicted instability is linear (exponentially growing), while for s = 2 it 
is nonlinear (sufexponential). 

From our numerical results, we have indeed found a weak instability mode with 
J = 1 for the s — 1 and s = 2 rings, see Fig. 5. To the limit of our numerical 
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Fig. 5. Weak instability of the ring solitons with s = 1 resulting from J = ±1 
perturbations. To the limit of the numerical method, this mode has zero instability 
growth rate at n > 0.16. A similar mode can be found for the s = 2 rings. 

accuracy, the corresponding instability growth rate ImA vanishes at k ~ 0.16, 
and then remains equal to zero up to the point k = ^offset ~ 0.18, at which the 
size of the ring becomes infinitely large, see Eq. (4). It is not completely clear 
yet whether an extremely small exponentially vanishing unstable eigenvalue 
exists past the value k ~ 0.16, but the issue is purely formal, the broad rings 
being stable in any practical sense. 

In conclusion, we have found that localised vortex solitons, with values of 
the spin 1 and 2, of the two-dimensional cubic-quintic nonlinear Schrodinger 
equation are linearly stable in a finite interval of the propagation constant 
corresponding to broad solitons with large power. These results make an im- 
portant step forward in the resolution of the controversial stability issue in 
this class of models. 
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